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Abstract 

This paper is concerned with the direct and inverse acoustic or electromagnetic scattering prob- 
lems by a locally perturbed, perfectly reflecting, infinite plane (which is called a locally rough surface 
in this paper). We propose a novel integral equation formulation for the direct scattering problem 
which is defined on a bounded curve (consisting of a bounded part of the infinite plane containing 
the local perturbation and the lower part of a circle) with two corners. This novel integral equation 
can be solved efficiently by using the Nystrom method with a graded mesh introduced previously by 
Kress and is capable of dealing with large wavenumber cases. For the inverse problem, we propose 
a Newton iteration method to reconstruct the local perturbation of the plane from multiple frequency 
far-field data, based on the novel integral equation formulation. Numerical examples are carried out 
to demonstrate that our reconstruction method is stable and accurate even for the case of multiple- 
scale profiles. 

Keywords: Integral equation, locally rough surface, inverse scattering problem, far field pattern, 
perfectly reflecting surface, Newton iteration. 

1 Introduction 

Consider problems of scattering of plane acoustic or electromagnetic waves by a locally perturbed, per- 
fectly reflecting, infinite plane (which is called a locally rough surface). Such problems occur in many 
applications such as radar, remote sensing, geophysics, medical imaging and nondestructive testing (see, 

e.g. m si 00 ed). 

In this paper we restrict the discussion to the two-dimensional case by assuming that the local per- 
turbation is invariant in the xj, direction. We assume throughout that the incident wave is time-harmonic 
{e~ l0Jt time dependence), so that the total wave field u satisfies the Helmholtz equation 

Au + k 2 u = in D+. (1.1) 

Here, D + := {{x\,X2) I *2 > hr(xi),x\ e R) represents a homogeneous medium above the locally rough 
surface denoted by Y := dD + = {(xi,X2) | X2 = hr(xi), x\ e R} with some smooth function hy € C 2 (R.) 
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having a compact support in R, k = w/c > is the wave number, w and c are the frequency and speed of 
the wave in D + , respectively. Throughout, we will assume that the incident field u' is the plane wave 

u\x) - exp(ikd ■ x), 

where d - (sin 9, - cos 9) € S _ is the incident direction, 9 is the angle of incidence, measured from the 
;t2-axis with -n/2 < 9 < n/2, and S _ := \x = {xi,x%) \ \x\ = \,x% <Q] is the lower part of the unit circle 
S = {x € R 2 | \x\ = 1}. We further assume that the total field u(x) = u'(x) + u'\x) + u s {x) vanishes on the 
surface T: 

u(x) = u\x) + u r (x) + u s {x) = on T, (1.2) 

where u r is the reflected wave by the infinite plane xj = 0: 

u r (x) = - exp(ik[xi sin 9 + X2 cos 9]) 

and u s is the unknown scattered wave to be determined which is required to satisfy the Sommerfeld 
radiation condition 

i (du s \ 

lrmrJ — - iku s = 0, r = \x\, x e D+. (1.3) 
i— x» \^ 5r / 

This problem models scattering of electromagnetic plane waves by a locally perturbed, perfectly conduct- 
ing, infinite plane in the TE polarization case; it also models acoustic scattering by a one-dimensional 
sound-soft surface. Figure [Dpresents the problem geometry. 

The well-posedness of the scattering problem (ll.lb - dl.3l l has been studied by using the variational 
method with a Dirichlet-to-Neumann (DtN) map in [4] or the integral equation method in QUI . In 
particular, it was proved in [30 1 that u s has the following asymptotic behavior at infinity: 

u\x) = — \u°°(x;d) + 0(-) 

uniformly for all observation directions x e 5+ with S + :- {x - (x\,X2) \ \x\ = l,X2 > 0} the upper 
part of the unit circle S , where u°°(x; d) is called the far field pattern of the scattered field u s , depending 
on the observation direction jc and the incident direction d e 5_. The integral equation formulation 
obtained in 11301 is of the second kind with a compact integral operator defined on the local perturbation 
part of the infinite plane. However, it is not suitable for numerical computation since it also involves an 
infinite integral over the unbounded, unperturbed part of the infinite plane. In [4], the scattering problem 
dl.lb - dl.3b is reformulated as an equivalent boundary value problem in a bounded domain with a DtN 
map on the part in D + of a large circle enclosing the local perturbation of the plane. This equivalent 
boundary value problem with a non-local boundary condition is then solved numerically by using the 
integral equation approach. However, the integral equation thus obtained involves a non-local DtN map 
on the semi-circle which needs to be truncated in numerical computations. 

In this paper, we propose a novel integral equation formulation for the scattering problem (|l.lM1.3b . 
which is defined on a bounded curve (consisting of a bounded part of the infinite plane containing the 
local perturbation and the lower part of a circle) with two corners. Compared with [4] and 11301 . our 
integral equation formulation does not involve any infinite integral or a DtN map and therefore leads 
to fast numerical solution of the scattering problem including the large wavenumber cases. In fact, our 
integral equation can be solved efficiently by using the Nystrom method with a graded mesh at the two 
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corners introduced previously by Kress [24] (see Section |3]below). Furthermore, we are also interested 
in the inverse problem of determining the locally rough surface from the far field pattern u°°(x, d) for all 
x e S + , d e S -. A Newton iteration method is presented to reconstruct the locally rough surface from 
multi-frequency far field data, and our novel integral equation is applied to solve the direct scattering 
problem in each iteration. From the numerical examples it is seen that multi-frequency data are necessary 
in order to get a stable and accurate reconstruction of the locally rough surface. 

The mathematical and computational aspects of the scattering problem (|l.U - (ll.3l have been studied 
extensively in the case when the local perturbation is below the infinite plane which is called the cavity 
problem (see, e.g. fl] |2] [6] and the references quoted there) and for the case of non-local perturbations 
which is called the rough surface scattering (see, e.g. ll8l l9l fTUl fTTl [T2l [T3l |3TT| >. 

There are many works concerning numerical solutions of the inverse problem of reconstructing the 
rough surfaces from the scattered field data. For example, a Newton method was proposed in ETl to 
reconstruct a local rough surface from the far-field pattern under the condition that the local perturbation 
is both star-like and over the infinite plane. An optimization method was introduced in [3] to recover a 
mild, local rough surface from the scattered field measured on a straight line within one wavelength above 
the local rough surface, under the assumption that the local perturbation is over the infinite plane. In (4), 
a continuation approach over the wave frequency was developed for reconstructing a general, local rough 
surface from the scattered field measured on an upper half-circle enclosing the local perturbation, based 
on the choice of the descent vector field. The reconstruction obtained in [4] is stable and accurate due to 
the use of multi-frequency near-field data (see also [5 1). It should be pointed out that the reconstruction 
algorithm developed in [4] does not work with multi-frequency far-field data. Note that our novel integral 
equation formulation can also be used to develop a similar Newton inversion algorithm with multiple 
frequency near-field data. For the numerical recovery of non-local rough surfaces we refer to []7] [14] [15] 
[19]|20]]. For the inverse cavity problem, the reader is referred to ll2l l2Tll28l . 

This paper is organized as follows. In Section [2] a novel integral equation formulation is proposed to 
solve the direct scattering problem. Section [3] is devoted to the numerical solution of the novel integral 
equation. In Section [4] it is proved that the local rough surface can be uniquely determined by the 
far-field pattern corresponding to a countably infinite number of incident plane waves. The Frechet 
differentiability is also shown of the far-field operator which maps the surface profile function hy to the 
corresponding far field pattern uT(x, d). The Newton method with multi-frequency far-field data is given 
in Section [5] based on the novel integral equation solver in Section [3] In Section [6] numerical examples 
are carried out to demonstrate that our reconstruction algorithm is stable and accurate even for the case of 
multiple-scale profiles, which is similar to the inversion algorithm with multi-frequency near-field data 
developed in H. 

2 A novel integral equation formulation for the direct problem 

Let / = -{u l + u r ). Then / is continuous on T and / = on To := {(jci, X2) £ T | X2 = 0}, that is, / has a 
compact support on T. The scattering problem (|l.lt - (ll.3t can be reformulated as the Dirichlet problem 
(DP) in the following way: 

Find u s € C 2 (D + ) n C(D+) satisfying the Helmholtz equation (11.11 ) in D + , the Sommerfeld radiation 
condition (11.31 ) and the Dirichlet boundary condition: 

u s — f on T. (2.1) 

The following uniqueness result has been proved in tl30l Theorem 3.1] for the above Dirichlet prob- 
lem (DP). 



3 




\3B 

i R r 



D 



Figure 1: The scattering problem from a locally rough surface 



Theorem 2.1. The problem (DP) has at most one solution in C 2 (D + ) n C(D + ). 

The existence of solutions to the problem (DP) has also been studied in [ 30 ] using an integral equa- 
tion method. However, the integral equation obtained in ll30l involves an infinite integral which yields 
difficulties in numerical computation. In this section, we propose a new integral equation to avoid this 
problem. To this end, we introduce the following notations. Let Br := {x = (x\,Xz) \ \x\ < R} be a circle 
with R > large enough so that the local perturbation T p := r\To = {(x\,hr(x\)) \ x\ e supp(/jr)l c Br. 
Then Tr := T D Br represents the part of T containing the local perturbation T p of the infinite plane. 
Denote by x A := (-R,0),x B :- (R,0) the endpoints of Tr. Write R+ := {(x\,xi) e R 2 | x 2 ^ 0}, 
D% := B R n D ± and dB% := dB R D ± , where D„ := {(xux 2 ) I x 2 < h T {x\), X\ e R}. See Figured] 

For (p € C(dD R ) define S k and to be the single- and double-layer potentials: 

(S k <p)(x) := f <D*(x,)0¥>O0<foG0. ieR 2 \^ 
(Djrtp)W := f 9 ^ ( , J ; y) y(y)^(y), xeR 2 \5D R 
and define S k,S r k e , K k , K[ e to be the boundary integral operators of the following form: 



(S k ifi)(x) := 



(srVX*) := 



I ® k (x,y)<p(y)ds(y), x e dD R 

JdD R 

I ® k (x,y)<p(y)ds(y), xeT R 
JdD- 



%(x re , y)0>OO<fo()O, xedB R U {x A , x B ) 



(K r k e cp)(x) 



■- f 

JdD- 



d® k (x,y) 
dv(y) 



(p(y)ds(y), x e <9D R 



1/ 

J 



d® k (x,y) 

bMy) 
d® k (x re ,y) 



<p(y)ds(y), xeT R 



dv(y) 



-tp(y)ds(y), x e dB R U {x A , x B } 

where x re = (x\,—X2) is the reflection of x = (x\,x 2 ) about the xi-axis, O k (x,y) is the fundamental 
solution of the Helmholtz equation Aw+k 2 w = with the wavenumber k, and v is the unit outward normal 
on dD R . Note that <t>o(x,y) = -1 /(2n) In \x - y\ is the fundamental solution of the Laplace equation. 
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Remark 2.2. Let <p e C(dD R ). From E2l Theorem 15.8b] it follows that the single-layer potential S k <p 
is continuous throughout R 2 . In addition, from ll26l Section 6.5] we know that the double-layer potential 
D (p can be continuously extended from R 2 \D R to R 2 \D R with the limiting value 



(K <pXx) + -<p{x) for x € 5D s \{x A , x b ) 

(K <p)(x) + — — (p{x) for x € {x A , x b ) 
2n 



where y(x) is the interior angle at the corner x e {xa, xb}- It remains valid for D k ip with k > since the 
kernel of D k - Do is weakly singular which yields that (D k - D$)<p is continuous throughout R 2 . Thus 
from the jump relations, S k , S r k e , K k and K r k are bounded in C(dD R ), where and K r , e are given by 

(K k (p)(x) for x € 5D~\{x A , x b ) 



:- i ( ^ >)(x) + (zg_l kW forx€{x A ,x B } 



(K k <p)(x) for x e <9D R \{x A , x B ) 

r(x) _ _n 

In 2 1 

In particular, S k - S o, S r k e - S r Q e ', K k - Ko and K r k - K™ are bounded in C(dD~). 

Let u s be the solution of the problem (DP). Then we can extend u s (x) into R 2 \Br by reflection, which 
is denoted again by u s {x), such that u s {x) = -u s {x re ) in R 2 \Br. By a regularity argument (see [16 page 
88] or (30] Theorem 3.1]) and the reflection principle, we know that u s e C 2 (R 2 \D R ) n C(R 2 \D R ) and 
satisfies the Helmholtz equation (11.11) in R 2 \D R . Following the idea in [ 17], we seek the solution u s in 
the form 

u s (x) = (Dk<p)(x) - ir](S k <p)(x), if e C(dD R ), x e R 2 \dD R (2.2) 
where r] + is a real coupling parameter. Let </r' e be a continuous mapping from dD R to <9D^ such that 



X, X € Tr 

x re , x e dfi^ U {xa, xb) 



Since w v (x) + u s (if/ re (x)) = -2[u'(x) + u r (x)] on Tr and w v (x) + u s (ifr re (x)) = on dB R U |x A , x#), and by 
the jump relations of the layer potentials, we obtain the boundary integral equation Ptp(x) = g(x), where 

( <p + (K k tp - ir]S k <p) + (K r k e <p - irjS r k e <p) , x e T R 

P<p := \ 1 , \ (2.3) 

^ -(p + (K k <p - irjS k <p) + \K[ e ip - irjS r k e (pj , x e dB R U |x A , x B ) 

and 

1 0, X€5B S U{X A ,X B } ^ ; 

Here, we have used the fact that the interior angles y(x) at the corners x A , xb are both n/2. Note that 
g e C(dD R ). Further, from Remark I2T21 and the continuity of \f/ re it follows that P is a bounded linear 
operator in C(dD R ). 

Conversely, we have the following result. 



5 



Lemma 2.3. Assume that u s is of the form ri2.2D with ip e C(dD R ) which satisfies the integral equation 
Pip = g with P, g defined in ( 12. 3\) , \2.4\ , respectively. Then u s € C 2 (D + ) n C(D + ) and solves the problem 
(DP). 

Proof. Since ip e C(dD R ), it follows from Remark |Z2]that u s e C 2 (R 2 \D R ) n C(R 2 \D~) and satisfies the 
Helmholtz equation dl.ll ) in R 2 \D S . In addition, Pep - g implies that u s (x) + u s (x re ) - on dB R U{x A , x B ] 
and u s (x) = -(u'(x) + u r (x)) on Y R . 

Let ii s (x) = -u s (x re ) in R 2 \B R . Then S v satisfies the Helmholtz equation dl.ll ) in R 2 \B R and S v (x) - 
m s (x) on <95#. Moreover, the uniqueness of the exterior Dirichlet problem (see, e.g. ifTTl Chapter 3]) 
implies that u s = u s in R 2 \B R . In particular, u s (x) - on T\Tr which yields that u s (x) = -{u l + u r ) on T. 
The proof is thus completed. □ 

We now prove the unique solvability of the integral equation Pip = g. 
Theorem 2.4. The integral equation Pip = g has a unique solution ip e C(dD R ) satisfying the estimate 

\\<P\\c(dD R ) < QW + u r \\ C(Y) (2.5) 

Proof. We need to deal with the corners. Following the idea in [ 171, we introduce the following boundary 
integral operators: for z = xa, xb 



Kg<p(x) 



f 

JdD R 

a 

Jdi 



dO (x,y) 
dv(y) 



tp(z)ds(y), x e dD R 



8D- 



f 

JdDZ 



dOo(x,y) 

dv(y) 
d®o(x re ,y) 

dv(y) 



ip(z)ds(y), xeT R 



ip(z)ds(y), x e dB R U {x A , x B ) 



For z - x A , xb define B £ (z) '■= (jc e R 2 | \x-z\ < e} with radius s small enough such that B e (z)n(r\r ) = 0. 
Choose a cut-off function^ € C^°(R 2 ) satisfying thatO <X < hx(x) = x(x re ),x = 1 in B e (xa) and^ = 
in B £ (xb). Since Ko tZ cp vanishes in R 2 \D R for z = xa, xb, we can rewrite (12.21 ) in the following form: 



u\x) = x(x) 



I ( 

JdDZ \ 



d® k (x,y) 
dv(y) 



ir]Q)k{x,y)\ip{y)ds{y) - J 

J JdDZ 



d®o(x,y) 
8v(y) 



ip(x A )ds(y) 



+ [1 ~X(x)] 



fl ( 

[JdD- \ 



dv(y) 
d®o(x,y) 



- iT)® k (x,y)\ip(y)ds(y) 



ip(x B )ds(y) 



x e 



-\D R 



Accordingly, using the jump relations of the layer potentials and the fact that x(x) = x(x re ), we rewrite 
Pip, defined in (I2.3I ). as Pip = I x ip + Aip + Bip, where 



(I x ip)(x) := 



ip(x), xeT R 
(l/2)[ip(x) +x<fi(x A ) + (1 -x)<P(x B )l xedB'U [x A , x B ] 



(Aip)(x) := -x(x)ip(x A ) - [1 -x(x)]<fi(x B ) 

(Bip)(x) := x(x) (K k ip - irjS k ip - K , XA ip) (x) + [1 ~x(x)] (K k ip - ii]S k ip - K ^ B ip) (x) 

+X(x) (K[ e ip - iqSfip - K r Q e XA ip) (x) + [1 - X (x)] (K£<p - iqS[ e ip - K^ip) (x) 
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Here, I x , A are bounded in C(dD R ). From Remark 12721 B is also bounded in C(dD R ). 
Step 1. We show that P is a Fredholm operator of index zero. 
Let 

(M <p)(x) := x(x) (K ip - K , XA <p) (x) + [1 ~x(x)] (K tp - K , XB (p) (x) 

+X (x) (K r e <p - K r Q * XA <p) (x) + [1 - X m (K r e <p - K r Q * XB <p) (x) 

From Remark I2T21 it is easy to see that Mq is bounded in C(dD R ). Since the integral operator B - Mo has 
a weakly singular kernel and \f/ re is a continuous mapping, the operator B - Mq is compact in C(dD R ). 
Thus, P - (I x + M ) = A + B - M is compact in C(dD R ). 

Moreover, for z = xa,xb and < r < e choose a cut-off function \j/ rA e (^(R 2 ) satisfying that 
< 1/^, < 1, ip r ,z( x ) = 1 in the region < \x - z\ < r/2 and tA r , z (x) = in the region r < \x - z\ < °°. 
Define M , r : C(dD R ) -» C(<9Dp by 

M , r <p := it/ r , XA x[K (il/ r , XA <p) - K , XA (if/ r , XA (p(x A ))\ 

+(Ar,x s (l -X^oOr,*^) - ^o.xsC'Ar,^^^))] 

+iff, %XA x[K r e (i// r ,x A <p) - K r e XA (^ XA ip(x A ))] 

+^, XB (l - x)[K r e (i/r r , XB <p) - K[( XB {^ XB ip{x B ))\ tp e C(dD R ). 

Since the kernel of Mo, r -Mo vanishes in a neighborhood of (x A , x A ) and (x^, jc^), it is compact in C(dD R ). 
Thus P - {I x + Mo >r ) is compact in C(dD R ) since P - (I x + Mo) is compact in C(dD R ). 
We now introduce the following norm on C(dD R ) : 

IML.O := max I max <p(*a))I + 1(1 -*)(</>- <£(*b)) I + \<f(x A )\ + \<p(x B )\] , 



max 

BB R U{x A ,x B } 



xX(<p-<p(.x A )) 



+ 



^(1 -x)(<P-<P(xb)) 



+ \<p{x A )\ + \tp(x B )\ 



which is equivalent to the maximum norm || • lU. It is easy to see that I x is a bijection from C(dD R ) to 
C(dD~) with 



- { 



1^, x e Fa 

2>p - x^(xa) - (1 - X)^{xb), x e d5 R U {x A , x B ) 



and ll/^^Hco < IML,o for all tp e C(dD R ). Furthermore, choose a function (f>o e C(dD R ) such that <po > 0, 
<Po(xa) = 4>o(x B ) = and <f>o reaches its maximum on Tr. Then H/^olU = H^olL.o, which implies that 

\\IxWc(dD R ,Ho,,o)^C(dD R ,\\-\U = !■ 

We now prove that l|Mo,rllc(50-,||-|U,o)^C(<?^,|HU) < 1 for r > small enough. For any ip € C(dD R ), 

supp(Mo ir </?) c B r (x A ) U B r (x B ), so we only need to consider (M^ r <p)(x) for x e dD R n (B r (x A ) U B r (x B )). 
Note first that for x e dD R n B r (^A) we have 

(M ,r<p)(x) = lf/ r , XA [K (l// r , XA <p) - K , XA (lff r ,x A <fi(x A ))] 

= if/ r , XA [K r Q e (^ XA (p) - K™ XA (lff r ,x A <f(X A ))] . 

Then the required estimate can be obtained by following the idea in ifTTl Section 3.5], together with the 
inequality: 

\v(y)-(x-y)\<C\x-y\ 2 (2.6) 
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for x e Tr U {xa, xg}, y £ T# or for x € <95#, y € <9B S . 
When x = xa, using (12.61 ) we have 

\(Mosp)(xa)\ < OlMko 

When x e Tr n fi r (xA)> the integral (Mo, r ^)(jc) can be split up into two parts: the first one over 
Tr n B,-{xa) and the second one over dB R n S,.(xa). The first part can be bounded with the upper 
bound Crll^Hoo o, estimated directly using (12.61 . Noting that v(y) • (x — y) does not change its sign for 
x e Tg n B r (jCA) and j € dB~ n Z?,-(xa), we have that for x € T# n B r (xA) the second part is bounded by 



JdB R nB r 



0*o(x,y) 



dv(y) 
d®o(x,y) 



(xa) 



2 
a(x) 



dB R nB r (x A ) 



dv(y) 

d®o(x,y) 
dv(y) 



tyr,x A ty){<p(y) - <p(x a )) 
^OOlMU.o 



ds(y) 



ds(y) 



IMU.o 



n 



IMIoo.O 



where a(x) is the angle between the two segments connecting x and the two endpoints of the arc dB R n 
B,-(xa)- Since the interior angle at x A is n/2, we have a(x) < 3n/4 if r is small enough. Thus, for 
x e Tr n B,-(xa) we have 

|(M ,^)(x)| < (Cr + 3/4)|ML, 

When x e dB~ n S,-(xaX the integral (Mo^Xx) can also be split up into two parts: the first one 
over Tr n B^(xa) and the second one over dB R n Z?,-(xa)- Since 5 In |x - y\/dv(y) = -X2/\x- y\ 2 for 
y € Tr n B,-(xa), then, for x € 5B R n B^Cxa) the first part is equal to 



t/fr,x A (x) | 
+<ffr >XA (x) I 



d®o(x,y) 



dv(y) r " XA 
d® Q (x re ,y) 
dv(y) 



^ r ,x A (y)(<p(y) - <p(x a )) 



ds(y) 



tr,x A (y)(<f(y) - <p(xa)) 



ds(y) = 0, 



whilst the second part can be estimated using (12.61 ) and is bounded with the upper bound OH^H^o- Thus, 
for x € 8Br^ n B,-(xa) we have 

|(M ,^)(x)| < O|ML t0 

Similar estimates can be obtained for (Mo ,-(fi)(x) when x e dD~ n 5 r (xg). Thus we have the estimate 

|(M () ,^)(x)| < (Cr + 3/4)IML. 

for x € dD~ n (B,-(xa) U B,-(xb)). Choosing r > small enough we obtain that 
ll^o,rllc(dD-,|HU,o)^C(dD-,|HU) < T Then, by the Neumann series, I x + Mo, r has a bounded inverse in 
C(dD~), yielding that P is a Fredholm operator of index zero since P - [P — (I x + Mo , )] + l x + Mo, r and 
P-(I X + M 0>r ) is compact in C(dD R ). 
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Step 2. We prove that P is injective and therefore invertible in C(dD R ). 

Let Ptp = for tp e C(dD R ). Then, by Lemma [2731 the scattered field u s defined by ( 12.21 ) vanishes 
on the boundary T. From Theorem 12. H and Holmgren's uniqueness theorem, u s vanishes in R. 2 \D^. On 
the other hand, by the jump relations of the layer potentials (see Remark |272l , we have (Pq<p)(x) = for 
x e dD R , where (Po<p)(x) '■- a(x)tp(x) + (K k tp)(x) - irj(Sk(p)(x) is bounded in C(dD R ) with a(x) = 1/2 
for x € dD R \{xA, xb) and a(x) = y(x)/(2n) for x = xa, *b- It was proved in [29] that Pq has a bounded 
inverse in C(dD R ), yielding that tp = 0. 

By the Fredholm alternative, P is invertible with a bounded inverse P~ l in C(dD R ). The proof is thus 
completed. □ 

Combining Lemma l2~31 Theorems 12. ll and[ 2.4l and the mapping properties of the single- and double- 
layer potentials in the space of continuous functions, we get the following result on the well-posedness 
of the problem (DP). 

Theorem 2.5. The problem (DP) has a unique solution u s € C 2 (D + ) n C(D + ). Furthermore, 



Remark 2.6. By the asymptotic behavior of the fundamental solution ®£, we have the following far field 
pattern for the scattered field u s given by (12.21) : 

lYl / a\ I rt y x A -i — lirY-V , x i , x it. 



which is an analytic function on the unit circle S , where <p is the solution of the integral equation Ptp - g. 

3 Numerical solution of the novel integral equation 

We make use of the Nystrom method with a graded mesh introduced in [ 17 , Section 3.5] (see also 11241 ) 
to solve the integral equation Ptp = g. Let dD R be parameterized as x(s) = (xi(s),X2(s)), < s < 2n 
such that (xi(0), *2(0)) = %b and (x\(jr), xzin)) = xa, where 



Here, on : [0, 2n] —> [0, 2n] is a strictly, monotonically increasing function satisfying that a>(s) - 
ji[v(s)] p / ([v(s)] p + \y(n - s)] p ) for0<s<n and a>(s) = a>(s - n) + n for n < s < 2n, where 



with p = 4. Note that a>(0) - 0, co(n) - n and oj'(0) = co'(7t) = 0. Then the integral equation P<p = g can 
be rewritten as Ptp{x(t)) = g(x(t)), where 



||« s || C(5+) < C||m ! + w r |lc(n- 



(2.7) 




(2.8) 





( r 



• d® k (x{t),x(s)) 

dv(x(s)) 
■ d® k (x(t),x(s)) 

dv(x(s)) 
d® k (x re (t),x(s)) 
dv(x(s)) 



(f(x(t))+2 I 



- i7]Q> k (x(t), x(s)) 



\x'(s)\ds, < t < n 



Jo 



P(p(x(t)) = \ -<p(x(t)) + 




- ir]<b k (x re (t), x(s)) \x(s)\ds, n<t<2n 



- ir]® k (x(t), x(s)) 



\x'(s)\ds 
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For j - 0, 1,. . . ,2n — 1 with n e M and n > 0, let - jn/n and = g{x{tj)). Then we get the 

approximation value of at the points x(t;) by solving the linear system 

j J 

-yf + 2 ^KMutfr + fatutfyf = g u i = 0,n 

2n-l 

<pf + 2 £ (R^K^j^ + ^itiJj^f = gi , i = l,..-, B -l 

2;i-l 



i = n + 1 , • • • , 2« - 1 



where, for j = 0, 1, . , . , 2n — 1 



1 n ^— j m \ n ) n 

m=\ 



mjn\ (-\) J n 



and for j = 0, 1, . . . ,2n — 1, j — 1, 2, • • • ,n — l,n + 1, • • • , 2n - 1 

"fl*jt(*(f«). *(*/)) 



Ki(ti,tj) 



dv(x(tj)) 



- ij]® k (x(ti), x(tj)) 



\x'( tj )\ 



1 

An 



dJ (k\x(td - X(tj)\) 



dv(x(tj)) 

K 2 (ti, tj) := K(t u tj) - K^U, tj) In (4 sin 2 
dO^itilxitj)) 



- iT]J (k\x(ti) - x(tj)\) 

tj - 1 



\ X '(tj)\ 



K 3 (ti, tj) := 



iT]® k (x re (ti), X(tj)) 



dv(x( tj )) 

Here, Jo is the Bessel function of order 0. Further, (I2.8t can be rewritten as 

a -inlA r 2n 



x\tj)\ 



u (x) = 



/Snk 



[k 

Jo 



v(x(s)) -x + tj] e- lkx - x(s) \x'{s)\ip(x(s))ds 



Then for rif £ IN with rif > we get the approximation values of the far field pattern at the points 
xi - in /rif, i = 0, 1, • • • , rif , by the following quadrature rule: 



„-ot/4 



2n-l 



M °°(x,)-^=-- V toi))-^ + ^"' Krl( ' J, k'(/,)l^ 



4 The inverse problem 

The inverse problem we are interested in is, given the far field pattern u°° of the scattered wave u s of the 
scattering problem fll.H - dl.3l (or the problem (DP)) for one or a finite number of incident plane waves 
u', to determine the unknown locally rough surface T (or the local perturbation T p ). 

We have the following uniqueness theorem which can be proved by arguing similarly as in the proof 
of Theorem 3.1 in ED- 
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Theorem 4.1. Assume that Y[ and Y2 are two locally rough surfaces and u°°(x, d) and u^{x, d) are the 
far field patterns corresponding to Y\ andY-i, respectively. Ifu™(x,d n ) - u^ '(x, d n ) for all x, € S + and 
d n e S- with n e N and a fixed wave number k, then Y\ - Yj- 

Given the incident plane wave u'(x) - exp(ikd • x) with the incident direction £ we define 
the far field operator F d mapping the function hr which describes the locally rough surface Y to the 
corresponding far field pattern u™(x,d) in L?(S +) of the scattered wave u s of the scattering problem 
(Ol-(fT31): 

F d {h T ) = u™{;d). (4.1) 

Here, we use the subscript k to indicate the dependence on the wave number k. In terms of this far 
field operator, given the far-field pattern u™(x,d), our inverse problem consists in solving the equation 
(14.11) for the unknown function hr- This is a nonlinear and very ill-posed operator equation. To solving 
this equation by the Newton method, we need the Frechet differentiability at hr- To this end, let Ah e 
C 2 R (R) := [h e C 2 (R) | supp(/i) c (-R,R)} be a small perturbation of the function h G € C 2 (R) and 
let r A /, := {(xi,hr(x\) + Ah(x\)) \ x\ € R} denote the corresponding boundary defined by /Jr(^i) + 
Ah(x\). Then F d is called Frechet differentiable at hy if there exists a linear bounded operator F'Jhr, •) : 
C 2 R (R) -> L 2 (S + ) such that 

\\F d (h r + Ah) - F d (h r ) - F' d {h T ,Ah)\\ LHs+) = o(\\Ah\\ C 2 m ) 

as ||A^|| C 2 (R) -> 0. 

Theorem 4.2. Let u(x,d) - u'(x,d) + u r (x,d) + u s (x,d), where u s solves the problem (DP) with the 
boundary data f - —{u l + u r ). If hi € C 2 , then Fj is Frechet differentiable at hy and the derivative 
F'AJit, Ah) = u'^for Ah e Cq R (R). Here, u'^ is the far field pattern of u' which solves the problem (DP) 
with the boundary data f = -(v2Ah)du/dv, where V2 is the second component of the unit normal v onY 
directed into the infinite domain D + . 

Proof. The proof is similar to that of Theorem 4.1 in [4 | with appropriate modifications. □ 

5 The Newton method with multi-frequency data 

We now describe the Newton iteration method for solving our inverse problem of reconstructing the 
function hi from the far field data, that is, for solving the equation (14. II ). Motivated by H, we use 
multi-frequency far field data in order to get an accurate reconstruction of the function hy. 

For each single frequency data with wave number k > 0, we replace (14.11 ) by the linearized equation 

F d (h r ) + F d (h T ,Ah) = u™(;d) (5.1) 

which we will solve for Ah by using the Levenberg-Marquardt algorithm (see, e.g. ll23l ) in order to 
improve an approximation to the function hi. The Newton method consists in iterating this procedure. 

In the numerical examples, we consider the noisy measurement data u™ k (x, di), x e S + ,l = l,...,n d , 
which satisfies 

\\uf k (;dd - U™{;di)\\ ms+) < 6\\u^{;dl)\y {s+ y 

Here, 8 is called a noisy ratio. In practical computations hi has to be taken from a finite-dimensional 
subspace Rm c: Cl „(R) and the equation (15.11) has to be approximately solved by projecting it on a 
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finite-dimensional subspace of L 2 (S +) by collocation at a finite number nf of equidistant points xj e S+, 
j = l,...,rif. Let Rm - span{(^i, 02, •• • ,4>m), where - 1,2, ... ,M, are spline functions with 
support in (-R,R) (see Remark I5TT1 below). Assume that h app e Rm is an approximation to hy with 
T app being the corresponding boundary. Then, by the strategy in ll23l . we seek an updated function 
Ah = Acii(f>j in Rm of h app such that Ac?,, / = 1, . . . , M, solve the minimization problem: 



(nd n f m 
Z \ F d,(h app )(xj) + F' d/ (h app , Ah){xj) - u^dtf +ySj] |AOi 
?=1 7 



(5.2) 



; = 1 



where the regularization parameter /3 > is chosen such that 



( n d n f y/ 2 

£ \F dl {h app ){ X j) + F' di {h app ,Ah){xj) - u^ k {xj,dt)\ 2 



P 



Y^lF^Xx^-u^dtf 



1/2 



(5.3) 



for a given constant p e (0, 1). Then a new approximation to /zp is given as h app + Ah. Define the error 
function 



2 n d 



Err k = — ) 



J |F rf ,(/^)(x ; ) - u^ddflj] K(*pdi)\ 2 



1/2 



Then the iteration is stopped if Err^ < t5, where t > 1 is a fixed constant. See ||23l for details. 

Remark 5.1. For a positive integer M € M + let /j = 2R/(M + 5) and = (/ + 2)/z - R. Then the spline 
basis functions of Rm are defined by 0,(?) = 0((£ - ti)/h), i — 1, 2, . . . , M, where 



k+i 



7=0 



c-iy 



+ 1 

j 



k+l 
t + ~ j 



with = z fc for z > and = for z < 0. In this paper, we choose k = 4, that is, is the cubic spline 
function. Note that 0,- e C 3 (R) with support in (-R,R). See HH for details. 

Remark 5.2. Our inversion Algorithm 15.11 below does not require the locally rough surface T to be 
parameterized by a function h r since, in practical computations, h^ is taken from the finite-dimensional 
subspace Rm spanned by spline functions (pi, j = 1,2, ...,M with support in {-R,R) (see discussions 
before Remark 1570 . Therefore, Algorithm 15.1 l ean deal with more general C 2 -smooth T. 

Remark 5.3. For the synthetic far-field data of the scattering problem, we choose the coupling parameter 
T] — k and get a finite number of measurements u^ k (xj,d), j = 0,1,..., rif, with equidistant points 
xj = jn/rif for a positive integer rif e N + . For the numerical solution of the scattering problem in 
each iteration, we choose t] - both to avoid the inverse crime and to reduce the complexity of the 
computation. Here, we need to assume that k is not a Dirichlet eigenvalue of the region bounded by 
the curves {(jci, h app {x{j) \ x\ € [-/?,/?]} and dB~. Further, it is seen from Theorem 14.2 1 that, in order 
to compute the Frechet derivative in each iteration we need to compute the normal derivative du s /dv of 
the scattered wave u s on a subset of {x{to), x(t\), . . . , x{t2 n -\)} which is contained in [{x\, h app {x\)) \ x\ € 
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(-R,R)}. Since h app € C 3 (R) and the two corners xa and xb are not included in the subset, we just use the 
quadrature rules in [25 1 in the form (12.21 of the scattered field u s and the graded mesh with the discrete 
values <p . , j = 0, 1, . . . , In - 1, of <p which are obtained from the scattering problem in each iteration. 

The Newton iteration algorithm with multi-frequency far-field data can be given in the following 
Algorithm 15. II 

Algorithm 5.1. Given the far field patterns u™ k {x,j,di),i = 1,2,. ..N,j - 0, 1, . ..,«/•, / = l,--- ,nj, 
where k\ < ki < ■ ■ ■ < k^. 

1 ) Let h app - be the initial guess ofh^ and set i — 0. 

2) Set i = i + \.Ifi>N, then stop the iteration; otherwise, set k - ki and go to Step 3). 

3) If Err \ < tS, return to Step 2); otherwise, go to Step 4). 

4) Solve (15. 2D with the strategy d5.il ) to get an updated function Ah. Let h app be updated by h app + Ah 

and go to Step 3 ). 

Remark 5.4. Our novel integral equation formulation proposed in Section [2] can also be applied to 
develop a similar Newton inversion algorithm with multiple frequency near-field data. 

6 Numerical examples 

In this section, several numerical experiments are presented to demonstrate the effectiveness of our algo- 
rithm. The following assumptions are made in all numerical experiments. 

1) For each example we use multi-frequency data with the wave numbers k = 1,3,..., 2N - 1, where 

Af is the total number of frequencies. 

2) To generate the synthetic data and to compute the Frechet derivative in each iteration, we solve the 

novel integral equation by choosing n — 128 for the wave number k < 13 and n = 256 for the wave 
number k > 13. 

3) We measure the half-aperture (the measurement angle is between and n) far-field pattern with 

65 measurement points, that is, nf = 64. The noisy data u™, are obtained as uT k = u™ + 
<%1|w£°|li2(s + )/llf Hi 2 (S + )' wnere f is a random number with Re(£), Im(£) € JV(0, 1). 

4) We set the parameters p = 0.8 and t - 1.5. 

5) In each figure, we use solid line '-' and dashed line '- -' to represent the actual curve and the recon- 

structed curve, respectively. 

6) For the shape of the local perturbation of the infinite plane in all examples, we assume that supp(/ir) € 

(-1, 1); we further choose R = 1 and use the smooth curves which are not in Rm- 

Example 1. In this example, we consider the case when the local perturbation of the infinite plane is 
over the x\ -axis with 

h T (xi) = (f>((xi + 0.2)/0.3), 
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Figure 2: The reconstructed curve (dashed line) at k = 1,5, 7, 11, respectively, from 3% noisy data with 
one incident direction d = (sin(7r/3), - cos(tt/3)), where the real curve is denoted by the solid line. 



where <p is defined in Remark 15.11 Here, we consider noisy data with 3% noise and use one incident 
direction d - (sin(7r/3), - cos(7r/3)). For the inverse problem, we choose the number of the spline basis 
functions to be M = lOand the total number of frequencies to be N = 6. Figure|2]shows the reconstructed 
curves at k = 1,5,7, 11, respectively. 

Example 2. In this example, we consider the case when the local perturbation of the infinite plane is 
under the xi-axis with 



where <p is also given in Remark [57TI In the inverse problem, the number of the spline basis functions 
is chosen to be M = 10 and the total number of frequencies is chosen to be N = 9. Figure [3] presents 
the reconstructed curves at k = 1, 5, 11, 17, respectively, from 3% noisy data with one incident direction 
d = (sin(7r/3), - cos(tt/3)). 

Example 3. The reconstruction considered in this example is a more challenging one with 



Here, we consider 10% noisy data. In order to get a good reconstruction, the number of the spline basis 
functions is taken to be M = 20 and the total number of frequencies is taken to be Af = 15. Figure [4] 
gives the reconstruction at k = 1, 9, 19, 29, respectively, with one incident direction d = (0,-1) (normal 
incidence from the top). 



h T (xi) = -O.80((xi - 0.3)/0.2), 
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k=11 





Figure 3: The reconstructed curve (dashed line) at k - 1, 5, 11, 17, respectively, from 3% noisy data with 
one incident direction d = (sin(7r/3), - cos(tt/3)), where the real curve is denoted by the solid line. 





-1 -0.5 0.5 1 -1 -0.5 0.5 1 



Figure 4: The reconstructed curve (dashed line) at k = 1,9, 19,29, respectively, from 10% noisy data 
with normal incidence from the top, where the real curve is denoted by the solid line. 
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Example 4 (multi-scale profile). We now consider the multi-scale case with 

( exp [l6/(25x? - 16)] [0.5 + 0.1 siii(l&rxi)] , |*i| < 4/5 
r(Xl) "I 0, |xi| > 0. 

This function has two scales: the macro-scale is represented by the function 0.5 exp[16/(25xj - 16)], 
and the micro-scale is represented by the function 0.1 exp[16/(25xj - 16)] sin(167rxi). To capture the 
two-scale features of the profile, the number of spline basis functions is chosen to be M = 40, and the 
total number of frequencies used is N = 30. The reconstruction is obtained with 10% noisy data using 
one incident plane wave with normal incidence from the top. Figure |5]presents the reconstructed profiles 
at k = 9,33,45,59. From Figure |5]it is observed that the macro-scale features are captured when k - 9 
(Figure [5] top left), while the micro-scale features are captured at k = 59 (Figure [5] bottom right). It is 
interesting to note that the resolution of the reconstruction does not improve much for k € [9, 33] and 
then improves from a larger k (e.g., k = 45) until a sufficiently large k (e.g., k = 59) for which the whole 
local rough surface is accurately recovered even with 10% noisy data. This indicates that our Newton 
algorithm with multiple frequency far-field data can give a stable and accurate reconstruction of multi- 
scale profiles with noise data as long as sufficiently high frequency data are used. This is similar to the 
reconstruction algorithm with multi-frequency near-field data developed in ||41 . 





-1 -0.5 0.5 1 -1 -0.5 0.5 1 



Figure 5: The reconstructed curve (dashed line) at k - 9,33,45,59, respectively, from 10% noisy data 
with normal incidence from the top, where the real curve is denoted by the solid line. 

The above numerical results illustrate that the Newton iteration algorithm with multiple frequency 
data gives a stable and accurate reconstruction of the local perturbation of the infinite plane even in 
the presence of 10% noise in measurements. From Figures EE it is seen that the upper part of the 
locally rough surface can be recovered easily at lower frequencies; however, much higher frequencies 
are needed in order to recover the deep, lower part of the locally rough surface as well as the fine details 
of the micro-scale features of multi-scale profiles. 
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We are currently trying to extend the technique to the TM polarization case. Furthermore, it is 
anticipated that the reconstruction method can be generalized to the three-dimensional case. 
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